#include "AnaDef.h"
#include "XGLUtils.h"
#include "dEdxGausN.h"

typedef Float_t vartype;
//typedef Double_t vartype;

vartype ftpcsig = -999;
vartype fpt = -999;

vartype feta = -999;  
vartype fv0cen = -999;

Float_t fdcaxy = -999;
Float_t fdcaz = -999;
Int_t ftpcncross = -999;

Bool_t DoCut()
{
  if(
     TMath::Abs(feta)<0.6 || TMath::Abs(feta)>0.8 || fv0cen>5 
     || TMath::Abs(fdcaz)>2 || TMath::Abs(fdcaxy)>0.5 
     || ftpcncross< 140 //10h, 10e cut
     ){
    return kFALSE;
  }
  else{
    return kTRUE;
  }
}

vartype HHx()
{ 
  //TPC ip large uncertainty in large p, not useful when converting to pt!! Use ESD p instead!!
  //return TMath::Log10(ftpcip); 
  //use ESD p

  //eta=-ln (tan (theta/2))
  //(1)
  //esd mom, pt, eta have
  //mom - pt/sin(2*atan(pow(TMath::E(),-eta)))
  //is always 0
  //(2)
  // for(int ii=0; ii<20; ii++){
  //   double pt= 10*gRandom->Rndm();
  //   double eta = -5+10*gRandom->Rndm(); 
  //   double esd = pt/sin(2*atan(pow(TMath::E(),-eta))); 
  //   TVector3 vec; vec.SetPtEtaPhi(pt, eta, 0); double my = vec.Mag(); 
  //   printf("%e %e --- %e \t\t %d\n", pt, eta, esd-my, TMath::Abs(esd-my)>1e-12);
  // }

  const vartype mom =  fpt/TMath::Sin(2*TMath::ATan(TMath::Power(TMath::E(),-feta)));
  return TMath::Log10(mom);
}

vartype HHy(){ return ftpcsig; }

vartype KKx(){ return HHx();}
vartype KKy(){ return TMath::Log10(fpt); }

void TestPrint(const Int_t id)
{
  printf("test%d hhx %e hhy %e -- %e %e %e %e %d\n", id, HHx(), HHy(), feta, fdcaxy, fdcaz, fv0cen, ftpcncross);
}

void extractHist(const TString fin)
{
  //input
  TChain *ch=XGLUtils::InputFiles(fin, "tree");

  ch->SetBranchStatus("*",0);

  ch->SetBranchStatus("tpcsig",1);
  ch->SetBranchStatus("pt",1);
  ch->SetBranchStatus("eta",1);
  ch->SetBranchStatus("v0cen",1);
  ch->SetBranchStatus("dcaz",1);
  ch->SetBranchStatus("dcaxy",1);
  ch->SetBranchStatus("tpcncross",1);

  ch->SetBranchAddress("tpcsig", &ftpcsig);
  ch->SetBranchAddress("pt",     &fpt);
  ch->SetBranchAddress("eta",     &feta);
  ch->SetBranchAddress("v0cen",     &fv0cen);
  ch->SetBranchAddress("dcaz",     &fdcaz);
  ch->SetBranchAddress("dcaxy",     &fdcaxy);
  ch->SetBranchAddress("tpcncross",     &ftpcncross);

  //ch->GetSelectedRows and ch->GetV1  do not work!!!!, many stange entries, maybe type conversion wrong!!

  TList *ll=new TList;

  //_____________________________________________________
  //_____________________________________________________

  TH2D * horiginal = dEdxGausN::GetHistTPCsigVsLog10P("horiginal");
  ll->Add(horiginal);

  TH2D * kk = dEdxGausN::GetHistLog10PtVsLog10P();
  ll->Add(kk);

  TH1D * heta = new TH1D("heta","",200,-1,1);
  ll->Add(heta);

  TH1D * hv0cen = new TH1D("hv0cen", "", 200, 0, 100);
  ll->Add(hv0cen);

  TH1D * hdcaxy = new TH1D("hdcaxy","",200, -2, 2);
  ll->Add(hdcaxy);

  TH1D *hdcaz = new TH1D("hdcaz", "", 200, -10,10);
  ll->Add(hdcaz);

  TH1I * htpcncross = new TH1I("htpcncross","",200,0,200);
  ll->Add(htpcncross);

  //original fill!
  Int_t nent = 0;
  while(ch->GetEntry(nent)){
    if(nent%1000000==0){
      printf("original fill %d\n", nent);
    }
    
    //test
    //TestPrint(0);

    if(DoCut()){
      horiginal->Fill(HHx(), HHy());
      kk->Fill(KKx(), KKy());

      heta->Fill(feta);
      hdcaxy->Fill(fdcaxy);
      hdcaz->Fill(fdcaz);
      hv0cen->Fill(fv0cen);
      htpcncross->Fill(ftpcncross);

      //test
      //TestPrint(1);exit(1);
  
    }

    nent++;
  }
  printf("original fill nent: %d\n", nent);

  //save original hyield
  TH1D *horiginalyields = dEdxGausN::GetHyield(horiginal);
  ll->Add(horiginalyields);
  
  //_____________________________________________________
  //_____________________________________________________

  ll->Print();
  ll->ls();

  TFile *fout=new TFile(Form("hist.root"), "recreate");
  ll->Write();
  fout->Save();
  fout->Close();

  printf("done!!");
}

int main (int argc, char * argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=2){
    printf("argc!=2\n"); return 1;
  }

  dEdxGausN::HistPrint();

  TStopwatch stwa;
  stwa.Start(kTRUE);

  extractHist(argv[1]);

  stwa.Stop(); 
  stwa.Print("m"); //"m": milli sec; "u": micro sec
  
  return 0;
}
